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The initial conditions describing seed fluctuations for the formation of structure in standard 
cosmological models, i.e. the Harrison-Zeldovich distribution, have very characteristic "super- 
homogeneous" properties: they are statistically translation invariant, isotropic, and the variance 
of the mass fluctuations in a region of volume V grows slower than V. We discuss the geometrical 
construction of distributions of points in R 3 with similar properties encountered in tiling and in 
statistical physics, e.g. the Gibbs distribution of a one-component system of charged particles in a 
uniform background (OCP). Modifications of the OCP can produce equilibrium correlations of the 
kind assumed in the cosmological context. We then describe how such systems can be used for the 
generation of initial conditions in gravitational iV-body simulations. 
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I. INTRODUCTION 

A central problem in contemporary cosmology is the 
quantitative explanation of the inhomogeneity observed 
in the Universe at large scales (see e.g. jl], ^|). Such in- 
homogeneity is probed both indirectly for early times via 
the fluctuations in temperature in the microwave sky, and 
directly at present in the distribution of matter in space. 
All currently standard cosmological models work within 
a paradigm in which these fluctuations are the result of 
(mostly) gravitational evolution operating on some initial 
very small fluctuations. Common to all such models is 
the assumption - supported in particular by observations 
of the microwave background Q- of a very specific form 
of these initial perturbations at large scale, i.e. Gaus- 
sian fluctuations with a Harrison-Zeldovich (HZ, often 
referred to as "scale- invariant" ) power spectrum |). 
In a recent paper some of us || have discussed how, in a 
simple classification of correlated systems, this spectrum 
corresponds to highly ordered glass-like or lattice- like dis- 
tributions. In particular such systems are characterised 
in real space by the striking behaviour of the variance in 
the mass, (AM 2 ), contained in a given volume, e.g. in 
spheres of radius R, (AM 2 ) oc R 2 , or i? 2 logi?, i.e. it is 
essentially proportional to the surface area of the sphere. 
This is to be contrasted with the Poisson behaviour (pro- 
portional to the volume R 3 ) characteristic of many equi- 
librium systems with short range interactions away from 
phase transitions or that characteristic of long-range cor- 
relations at a critical point in many systems (for which 
(AM 2 ) oc R a with a > 3). Actually R 2 is the slowest 



growth possible for tmwisotropic translationally invariant 
distribution of points 0] ■ 

In this paper we examine the nature of the correlations 
in such super-homogeneous systems and compare them 
with equilibrium correlations in systems with long range 
interactions studied in statistical mechanics p3[ . 
More precisely we consider the One Component Plasma 
(OCP), which is simply a system of charged point par- 
ticles interacting through a repulsive 1/r potential, in a 
uniform background which gives overall charge neutrality 
(for a review, see JllJ ) . We discuss how the OCP can be 
modified so that its equilibrium correlations are precisely 
those considered in cosmology. This involves an appro- 
priately designed attractive short range potential as well 
as a repulsive 1/r 2 potential at large scales. The lat- 
ter corresponds to a four dimensional Coulomb potential 
with the particles confined to three dimensions. 

Our analysis of these parallels has more than theoret- 
ical interest. It also provides a new method for gener- 
ating initial conditions (IC) for numerical studies of the 
formation of structure in cosmology [0, [l3| . We explain 
that this method avoids the problems associated with 
currently used algorithms which involve superimposing 
small perturbations on a lattice or glass-like distribution. 
The latter is understood to be a "sufficiently uniform" 
discretization of a constant density background, not as 
representing the perturbed distribution. The problem 
with this choice of initial state is that the unperturbed dis- 
tribution has its own inherent fluctuations/correlations, 
which can play an important role in the evolution of the 
system @ [|,[f],0. 
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II. THE HZ SPECTRUM IN COSMOLOGY 

Let us recall first the necessary essentials of modeling 
perturbations from a uniform density used in cosmology. 
These are described as a continuous stationary stochastic 
process (SSP), with the property that the volume aver- 
age of observables, i.e. empirical averages, approach the 
ensemble averages in the large volume limit. The sta- 
tionarity refers in this context to the statistical invari- 
ance under spatial translation. Moreover one assumes 
also statistical isotropy of the stochastic process, i.e. in- 
variance under rotation of the ensemble quantities. Thus 
if p(r) is the density field, we have 



(p(r)) = Po > , 



(1) 



and we can define a reduced 2-point correlation function 
fa) by 



(p(r\)p(r 2 )) 



Pi 



(2) 



where r = \r\ — Corresponding isotropies are as- 
sumed to hold for the higher order correlation functions. 
Further the probability distribution for the fluctuations 
in the initial uniform density are generally assumed to be 
Gaussian. 

Much more used in cosmology than £(r) is the equiva- 
lent fc-space quantity, the power spectrum P(k) which is 
defined as 



P(k) 



lim 



\Sp(k)\ : 



(3) 



where S p (k) = f v d 3 re~ lk r 5(r) is the Fourier integral 
in the volume V of the normalized fluctuation field 
5(f) — (p(r) — po)/pa. In a statistically isotropic SSP 
this depends only on |fc| with 



P{k) = P(k) 



£(r) — p — '-4irr 2 dr. (4) 
kr 



It follows from its definition that P{k) > 0, and, in the 
case we are considering, its integral, which is equal to 
£(0), is assumed to be finite. 

In current cosmological models it is generally assumed 
that the power spectrum P(k) cx fc at small fc. This 
is known as the HZ or "scale-invariant" power spec- 
trum. It is believed to describe the "primordial" fluc- 
tuations at very early times, the putative remnants of 
a period of "inflation" prior to the ordinary Big Bang 
phase (ij D . The linear fc behaviour is cut-off by a short 
distance scale (needed to ensure integrability) i.e. for fc 
larger than some cut-off, k c , P(k) decreases faster than 
fc -3 . The reason for the appearance of this particular 
spectrum is tied to considerations about the cosmolog- 
ical model (see ||). Thus in the homogeneous and 
isotropic Friedmann-Robertson- Walker (FRW) cosmol- 
ogy, the only characteristic length scale is the size of the 



horizon Rn(t) (the region causally connected at a time t). 
The HZ spectrum then corresponds to the choice which 
gives (as discussed in M this condition is in fact only sat- 
isfied for a spherical Gaussian window function, and not 
in a real physical sphere) 



Oj4 {Rh (t)) — constant, 



(5) 



independent of t, where a\ I (R) is the normalized mass 



variance 
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(R) = (AM 2 ) I (MY 



(6) 



The HZ spectrum then follows from a consistency crite- 
rion for the treatment of the perturbations in the model: 
any other spectrum will give mass fluctuations which 
dominate over the homogeneous background either at 
some time in the future or past. 

All current cosmological models share the HZ spectrum 
(or something very close to it) as the "primordial" form of 
their perturbations. This spectrum evolves, at larger and 
larger scales, until the characteristic time at which the 
densities of matter and radiation energy are equal. This 
evolution which depends strongly on the particular model 
used may change completely the form of the spectrum 
for distances smaller than a characteristic scale k~ q x (i.e. 
k > k eq ). For larger physical scales (k < k eq ), however, 
the evolution is almost exactly the same in all models, 
and it leaves the primordial HZ form intact. 

Viewed in the general framework of correlated pro- 
cesses a very crucial property of the HZ type power spec- 
trum is that 



lim P(k) = 
which in real space implies that 



(7) 



lim - — — — - = lim 

R^oc V(R) R^oo 



d 3 r i(r) = (8) 



V(R) 



where V(R) is the volume of a sphere of radius R (with 
arbitrary origin). Eq. (g) is to be contrasted with the 
behaviour of this integral in a Poisson distribution, when 
it yields a finite positive constant, and with that asso- 
ciated with power-law correlations in critical systems || 
for which the same integral of £(r) diverges. 

Systems satisfying Eq. (||) are thus "more homo- 
geneous" than a Poisson type system, as can be seen 
when one considers the behaviour of the mass variance 
in spheres. For a spectrum such that P(k) ~ fc™ for 
k — ► and appropriately cut-off at large k we have for 
large R 



' M 



(i?) (X 



l/i? 3+n if n< 1 
log(ii) /R 4 if n = 1 
l/R 4 if n > 1 



(9) 



In terms of the non-normalized quantity ((AM) 2 ) cx 
a 2 j(R)R 6 , Eq. (8) says that for n > we have a slower 



FIG. 2: Fragmentation step for the "pinwheel" tiling of the 
plane.) 



FIG. 1: A super-homogeneous distribution (bottom) and a 
Poisson distribution (top) with approximately the same num- 
ber of points. Both are projections of thin slices of three 
dimensional distributions. The super-homogeneous one is a 
"pre-initial" configuration used in setting up the configura- 
tion representing the initial conditions for the cosmological 
N-body simulations of j^]. As described in section [vi] the dy- 
namics used to generate this "pre-initial" configuration gives 
rise to surface fluctuations as it is essentially the one compo- 
nent plasma which we discuss below. 

increase of the variance as a function of R than for Pois- 
son fluctuations, corresponding to n = 0, with the limit- 
ing behaviour corresponding to a variance which is pro- 
portional to the surface area of the sphere. These sys- 
tems are thus characterised by surface fluctuations, or- 
dered (or homogeneous enough, one could say) to give 
this very particular behaviour. In Fig. 1 are shown, in 
two dimensions, a Poissonian distribution and a super- 
homogeneous distribution. The relatively greater unifor- 
mity of the latter is clearly identifiable. 



that in the case of a lattice is the square root of the 
average number of points in the last spherical shell of 
thickness equal to the lattice unit. A better (more sta- 
tistically uniform) example of the same kind is the so- 
called shuffled lattice ||; this is a lattice whose sites are 
independently randomly displaced by a distance x in all 
directions from their initial position according to some 
distribution p{x) which has a finite second moment. In 
this case we find P(k) ~ k 2 at small k and, consequently, 
again a 2 A {R) cx 1/i? 4 at large R. The simple lattice, how- 
ever, is not a SSP, and even the shuffled lattice, though it 
can be defined as a SSP, is not statistically isotropic be- 
cause the underlying lattice structure is not completely 
erased by the shuffling || . 

To construct a statistically isotropic and homogeneous 
particle distribution with such a behaviour of <J 2 {R) is 
non-trivial. A particular example is the so-called "pin- 
wheel" tiling of the plane |l9|, [2(J. The generation al- 
gorithm for it is defined by taking a right angled tri- 
angle with sides of respective length one and two (and 
hypotenuse a/5) and, at the first step, forming five sim- 
ilar square triangle of sides 1 / a/5 and 2 / a/5 respectively 
as shown in Fig. [|. At the second step we expand these 
new triangles, to the size of the original triangle, and 
repeat the procedure ad infinitum, so that they cover 
the plane completely. Finally, placing a point randomly 
inside each elementary triangle will give the superhomo- 
geneous point distribution which is statistically isotropic 
(with a continuous power spectrum). 



III. DISTRIBUTIONS OF POINTS WITH 
SURFACE GROWTH OF THE SQUARED 
VARIANCE 

A simple cubic lattice is the simplest example of a dis- 
crete set of points which shows this limiting behaviour 
(aftf oc 1/i? 4 ) of the variance |l^]. The result is not hard 
to understand: The standard deviation ctm is that mea- 
sured by averaging over spheres of radius R centered at 
a randomly chosen point of the unit cell. It is propor- 
tional to the typical variation of points in the volume, 



IV. A PHYSICAL EXAMPLE: THE 
ONE-COMPONENT PLASMA 

In the study of the one component plasma (OCP) one 
considers the equilibrium statistical mechanics (at inverse 
temperature 0) of a system of charged point particles 
interacting through a Coulomb potential, in a continuous 
uniform background giving overall charge neutrality JTl]| . 
Taking the particles to carry unit mass and charge the 
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microscopic mass density of the particles is given by 

p(f)=^2S(P-fi) (10) 

i 

where is the position of the ith particle. The total 
microscopic charge density is 



(ii) 



where 

n = (p(r)) . (12) 
The two-point correlation function is defined as 



(p) 2 

6(f-f") 



+ h(r- f) . 



(13) 



The first term comes from the diagonal terms r = r ' in 
(p(r)p(r')) and is characteristic of any point distribution. 

It can be proven quite generally H that at high tem- 
peratures the OCP has a translation invariant isotropic 
distribution with £(r) satisfying Eq. (||) with r = \r\. In 
terms of h(r) this gives 



if h(r)d 3 r = -1 



(14) 



Eq.(pT[) is usually referred to as a sum rule which origi- 
nates from the perfect screening of each charge caused by 
the long-range nature of the Coulomb potential. To see 
this let us suppose that an external infinitesimal charge 
density pext = ee lfc r of very long wavelength is applied 
to the system. It creates an external electric potential 



4tt 



1 ' J \f-f'\ * 2 
and a perturbation to the Hamiltonian 



(15) 



H cxt = J p c (f)$(f)d 3 r= f p c (r)e^ p d 3 r. (16) 

Using linear response theory the induced charge in the 
system is given by 



Pind(r') = ~(3(p c {f')H ext ) 



(17) 



where the average is over the unperturbed statistical dis- 
tribution. Thus, assuming that the applied charge is per- 
fectly screened i.e. the system responds with an induced 
charge density pi n( j = —p C xt we have, in the limit, k — > 0, 



ik.r' 



-Itt 



fcjj / ( Pc (r')Pc(?)W k - r d 3 r (18) 



and therefore, since 



where h(k) is the Fourier transform of h, we find that 



P(k) 



Aim 2 13 



(20) 



for small k. 

The behaviour of the power spectrum at small k is 
traceable through Eq. ( |l5| ) as being simply, up to a fac- 
tor /3, the inverse of the Fourier transform, (in the sense of 
distributions) of the repulsive 1/r potential. It is evident 
that what one needs to obtain the same kind of behaviour 
but with P{k) ~ k at small k is a repulsive potential of 
which the Fourier transform behaves as 1/fc, i.e. a re- 
pulsive 1/r 2 potential. Before considering this in more 
detail we discuss some issues relevant to the cosmological 
context which we have passed over without comment so 
far. 



V. DISCRETE VS. CONTINUOUS 

An essential difference between the case of cosmolog- 
ical perturbations and the system just discussed is that 
the former refers to a continuous density field while the 
latter describes the correlation properties of a set of dis- 
crete points. To interpret the latter as giving a realiza- 
tion of the former we obviously need to specify explicitly 
how to relate the two. There is no unique prescription to 
pass from a discrete to a continuous field. A prescription 
simply corresponds to a regularisation of the Dirac delta 
with a function VFi(r) with the property 



W L {r) = L~ 3 W (-), 



W {r)d 3 f= 1 (21) 



where L is the characteristic scale introduced by the reg- 
ularisation e.g. the Gaussian 



W L (r) 



1 



2nL 



exp 



r 

2I 2 " 



(22) 



For any finite value of L we can define a continuous den- 
sity field pl(t) as the convolution of this function with 
the discrete density field. The pair correlation function of 
the continuous field can then be written also as a double 
convolution integral of Wi,{r) and the correlation func- 
tion in the discrete case. The singularity in the latter 
is thus also removed by the regularisation. The power 
spectrum of this continuous field is then simply given as 



P L (k) = \W L (k)\ 2 P D (k) 



(23) 



where WL(k) is the Fourier transform of the regulari- 
sation Wi(f) and Po(fc) is the power spectrum of the 
discrete field. In particular for the Gaussian smoothing 
function wc have 



P L (k) =exp(-k*L J )P D (k). 



(24) 



P(k) = -+h(k) 
n 



(19) 
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For the OCP discussed above the correlation function, 
in the weak coupling limit, is approximately the Debye- 
Hiickel formula 



h(r) = -fl- 



it = 4n(3n 



so that 



P D (k) = 



K 2 + k 2 



(25) 



(26) 



which is a non-integrable function of k. Integrability is 
ensured by the cut-off imposed at large k in the power 
spectrum through the smoothing function. An assumed 
property of the primordial density field in standard cos- 
mological models is that it is Gaussian distributed. Since 
the density field is inherently positive this assumption of 
Gaussianity is more properly attributed to fluctuations 
(expected to be small) around the uniform density. Small 
fluctuations in the discrete OCP are in fact Gaussian [E2|. 



VI. INITIAL CONDITIONS IN N-BODY 
SIMULATIONS 

A context in which it is necessary to have a concrete 
realization of the density p(r) in cosmological models is 
that of the N-body approach to the problem of struc- 
ture formation ]l2] , In these numerical simulations 
an initial configuration (1C), which should represent the 
universe at early times, is evolved under its self-gravity. 
The goal of these simulations is an understanding of how 
structures grow in such models, and whether they are 
compatible with observations of the distribution of mat- 
ter (as probed, primarily, by the distribution of galaxies 
and galaxy clusters) . The initial continuous density field 
is thus necessarily represented by a discrete set of points, 
which, because the discretization scale does not represent 
the real underlying physical one, must be considered just 
as we have discussed above as representing the continu- 
ous density field via some appropriate smoothing. 

In practice such a discretization is always generated in 
a very particular way Q : points of equal mass are 
displaced from a "pre-initial" configuration in a way pre- 
scribed by the correlation field one wants to set up. This 
prescription can be understood as follows. Superimpos- 
ing an infinitesimal displacement field u(r) on a uniform 
density p Q one has 



P(r) - Po 



Po 



-V.u(f) 



(27) 



and therefore in /c-space the power spectrum of the cor- 
related density field p(f) produced is roughly propor- 
tional to k 2 P u (k), where P u (k) is that of the displace- 
ment field. The "pre-initial" configuration used, which 
gives the initial unperturbed positions of the particles, 
is a discretization of the initial uniform density field p 
and the displacement field is specified by the power spec- 
trum of the continuous field which is desired as IC in the 



simulation. In practice the former is taken either to be 
an exact simple lattice or a "glassy" configuration |fl3| . 
When this latter configuration is used it is generated by 
evolving the system with the N-body code, but with the 
sign of gravity reversed. This corresponds to a system 
essentially like the OCP (at low temperature), and the 
long-time evolution brings it to a highly uniform configu- 
ration with normalized variance decaying at large scales 
asf/i? 4 Jl|. 

In the generation of the IC in this way the fluctua- 
tions associated with the "pre-initial" configuration are 
neglected, or rather implicitly assumed to be negligible 
with respect to those which are introduced by the per- 
turbations to it. The problem with this procedure is 
that it neglects precisely these "pre-initial" fluctuations. 
(If we assume the initial continuous field to have power- 
spectrum Pi(k) what results from the infinitesimal shuf- 
fling is a power spectrum P(k) = Pi(k) + k 2 P u (k)). Such 
a neglect is justified if these fluctuations play no role in 
the evolution of the system, as one wants to see only the 
growth of structure coming from the perturbations su- 
perimposed on the uniform density field. Indeed the rea- 
son the lattice or "glassy" configuration is used is that 
they are configurations in which the gravitational force 
is (effectively) zero, whereas in a Poisson distribution for 
example (which is a priori as good a discretization of the 
uniform background) it is not. While the Poisson distri- 
bution will thus evolve and form structures even without 
additional perturbation, the perfect lattice will not. That 
the discretized lattice or glass solves this problem is how- 
ever far from clear: they both represent unstable point 
configurations with respect to gravity. A small applied 
perturbation will thus in general give an evolution which 
depends on the correlations and/or fluctuations in the 
initial unperturbed configurations. In principle such ef- 
fects can be kept under control by making the discretiza- 
tion of the original density sufficiently fine. In practice 
in N-body simulations the "non-linear regime" in which 
structure formation takes place is comparable to the dis- 
cretization scale, and the role of such effects seems, at 
the very least, problematic jjj, U, 0- 

One way in which this difficulty could be got around 
is by generating directly a discrete distribution whose 
regularisation is exactly the desired continuous field. Our 
observations on the OCP provides in principle a way of 
realizing this possibility. 



VII. IC AS EQUILIBRIUM OF A MODIFIED 
OCP 



We have seen that the OCP equilibrium correlations 
give surface fluctuations ((AM 2 ) ~ R 2 ), but with a 
power spectrum at small k which goes like k 2 . By con- 
sidering instead a repulsive 1/r 2 potential, whose Fourier 
transform is 2n 2 /k, we obtain from Eq. ( p9| ) 



P(k) 



2n 2 n 2 P 



, for k ~ 



(28) 
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and 



h(r) ~ - 



27r 4 /3n 2 



(29) 



for r 3> (27r 2 /3m) _1 . The change from the exponential 
decay of Eq. (|25| ) to a power-law decay is a result of the 
different analyticity properties of the two power-spectra: 
the k 2 behaviour is analytic at the origin, guaranteeing 
a rapidly decaying behaviour of its Fourier transform, 
while the k spectrum is not. 

In the context of cosmological N-body simulations 
what one needs is not simply the primordial HZ power 
spectrum with some appropriate small scale cut-off: what 
is simulated is only a part of the cosmological evolu- 
tion, starting from a time at which the initial spectrum 
of fluctuations is already significantly modified from its 
primordial form. While purely gravitational evolution at 
these early times does not modify the HZ spectrum, non- 
gravitational effects, present until the time when the uni- 
verse becomes dominated by matter, do so significantly. 
The nature of these modifications depends on the de- 
tails of the cosmological model, but in all cases it affects 
only k larger than a characteristic k eq corresponding to 
the causal horizon at the "time of equality" (when mat- 
ter starts to dominate over radiation). One then has a 
power spectrum of the form 



P(k) = k n f(k) 



(30) 



with n exactly (or very close to) unity and f(k) a model- 
dependent form for the spectrum with /(0) = a > 0. In 
Figure |^ are shown the power spectra of two represen- 
tative cosmological models |W ^[ . Both are CDM ( "cold 
dark matter") type. HDM ^'hot dark matter") models, 
now strongly disfavored, have a more abrupt (typically 
exponential) cut-off at larger k, due to the "wash-out" of 
small scale structure associated with their higher velocity 
dispersions. 

To produce such a spectrum as the equilibrium one of 
an OCP like system requires further modification of the 
form of the potential. Just like the standard OCP, an 
unmodified 1/r 2 potential will give, in the weak-coupling 
limit, a spectrum which becomes flat (i.e. Poisson like) 
at large k. This is simply due to the fact that in this 
limit the thermal fluctuations dominate over the poten- 
tial at small scales, effectively randomizing the particle 
positions up to some scale. A crude guess of what poten- 
tial would produce the behaviour of a typical cosmolog- 
ical model can be obtained by supposing that at small 
scales 



I + h(r) « e 



-0V{r) 



(31) 



which corresponds to completely neglecting collective ef- 
fects (1 + h represents the relative probability compared 
to random of finding a particle at distance r from a 
given one). Given that the desired fluctuations always 
have small amplitudes (\h\ <C 1), we would then need 
to be in the regime of temperature and scales such that 
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FIG. 3: Power spectra of two different cosmological models, 
both of the CDM ("cold dark matter") type but with a dif- 
fering density of matter (parametrised by Q). The plot is in 
log-log units so that the behaviour P(k) ~ k and P(k) ~ k~ 2 
correspond to straight lines of slope +1 and —2 respectively. 



\PV(r)\ « 1, so that h(r) /3V(r). The potential 

should thus be attractive at smaller scales, as the system 
is positively correlated at those scales. A k~ 2 behaviour 
at large k, which is often used (see Figure ||) as an ini- 
tial condition approximating cosmological models (CDM 
type) in this regime (beyond the "turn-over"), may be 
obtained from an attractive 1/r potential, e.g. 



(32) 



By modifying the parameters \x and rj, as well as the 
temperature, both the amplitude of the P(k) and the lo- 
cation of a change from a P(k) ~ k~ 2 to a behaviour 
P(k) ~ k for small k can be controlled. The potential 
( |32| ) is repulsive at short distances, but it may be nec- 
essary to make it more strongly repulsive in order to to 
ensure the system is not unstable to collapse [25 . Such 



issues, as well as the practical numerical implementation 
of these observations will be investigated in forthcoming 
work. Once this method has been implemented to set 
up the initial conditions, it will naturally be interesting 
to understand the effect of this change in initial condi- 
tions on the dynamical evolution of gravitational N-body 
simulations. 
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